#######################################################

library(pheatmap)
library(data.table)
library(optparse)

##########################################################################################
option_list <- list(
    make_option(c("--exp_file"), type = "character"),
    make_option(c("--diff_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 输出
    out_path <- "~/20231121_singleMuti/results/celltype_plot/diff_expression/germ/plot/"

    ## 表达文件
    exp_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ/GeneExpression.MeanByCellType.tsv"

    ## 
    diff_file <-"~/20231121_singleMuti/results/celltype_plot/diff_expression/germ/one_vs_other.germ.pct_0.25.tsv"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

diff_file <- opt$diff_file
exp_file <- opt$exp_file
out_path <- opt$out_path

dir.create( out_path , recursive = T )

###########################################################################################
## 读数据
a <- read.table(diff_file,sep="\t",header=TRUE)
dat_exp <- data.frame(fread(exp_file))

## 提取差异表达基因
a_subset<-subset(a,a$avg_log2FC>1&a$p_val_adj < 0.05)
rownames(dat_exp) <- dat_exp$gene

###########################################################################################
## 细胞顺序
grp_order2 = c("SSC",
"Differenting&Differented SPG",
"Leptotene",
"Zygotene",
"Patchytene",
"Diplotene",
"Early stage of spermatids",
"Round&ElongateS.tids",
"Sperm",
"Leydig cells",
"Myoid cells",
"Pericytes",
"Sertoli cells",
"Endothelial cells",
"NKT cells",
"Macrophages")

cell_list_germ_combine <- c("SSC_SPG" , "SPC" , "SPT")

cell_list <- c(grp_order2 , cell_list_germ_combine)

## 提取用到的细胞类型
cell_list <- cell_list[cell_list %in% unique(a_subset$cell_compare) ]

###########################################################################################
## 按照foldchange排序
all_gene_list<-c()

for(cellN in cell_list){

  print(cellN)
  data_cell <- subset( a_subset , cell_compare == cellN )
  sorted_data_cell <- data_cell[order(data_cell$avg_log2FC, decreasing = TRUE), ]
  gene_list <- sorted_data_cell$gene
  all_gene_list<-c(all_gene_list,gene_list)
}

un_all_gene_list<-unique(all_gene_list)

###########################################################################################
## 保留差异基因,修改细胞顺序
if( length(grep("&" , cell_list)) > 0 ){
  cell_list <- gsub( "&" , "." , gsub( " " , "." , cell_list) )
}else if( length(grep(" " , cell_list)) > 0 ){
  cell_list <- gsub( " " , "." , cell_list)
}

b_subset <- dat_exp[un_all_gene_list,cell_list]

p <- pheatmap(b_subset,
  scale = "row",
  show_rownames = F ,show_colnames = T, 
  #cutree_cols=7,cutree_rows = 7,
  cluster_rows = F , cluster_cols = F , clustering_method = "ward.D2",
  color = colorRampPalette(c("blue", "white", "firebrick3"))(100)
  )

out_file <- paste0( out_path , "/" , "heatmap.pdf" )
pdf(out_file , width = 8, height = 12)
print(p)
dev.off()
